#include "KSimPrimaryGeneratorAction.hh"
#include "KSimGlobals.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"
#include <G4ParticleTable.hh>
#include <G4ParticleDefinition.hh>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cmath>
using namespace std;

KSimPrimaryGeneratorAction::KSimPrimaryGeneratorAction(int neutronCount) {
	particleGun = new G4ParticleGun(1);
	G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
	G4ParticleDefinition *particle;
	particle = particleTable->FindParticle("neutron");
	particleGun->SetParticleDefinition(particle);
	particleGun->SetParticleTime(0);
	numOfNeutrons = neutronCount;
	randEngine = new CLHEP::MTwistEngine(48579284);
	cout<<endl;

}

KSimPrimaryGeneratorAction::~KSimPrimaryGeneratorAction() {
	delete particleGun;
	delete randEngine;
}

void KSimPrimaryGeneratorAction::GeneratePrimaries(G4Event* event) {
	eventNum++;

	double z = -6.15 * m;
	particleGun->SetParticlePosition(G4ThreeVector(0,0,z));
	double polLow = atan(.141421356 * m / (-z));
	double polarAngle = (atan(.9 * m / (-z)) - polLow) * randEngine->flat() + polLow;
	double azimuthalAngle = randEngine->flat() * 6.28318531;
	particleGun->SetParticleMomentumDirection(G4ThreeVector(cos(azimuthalAngle) * sin(polarAngle), sin(azimuthalAngle) * sin(polarAngle), cos(polarAngle)));
	neutronEnergy = (randEngine->flat()) * MeV;
	particleGun->SetParticleEnergy(neutronEnergy);
	particleGun->GeneratePrimaryVertex(event);

	stringstream ss;
	ss.width(15);
	ss.fill('0');
	ss << left << ((double) eventNum) / ((double) numOfNeutrons) * 100.0;
	cout << "\r" << ss.str() << " % done" << flush;
}
